#include "topology.hpp"
#include "meshTypes.hpp"
#include "signalCodes.hpp"
#include "util.hpp"
#include "geometry.hpp"
#include "streamAccess.hpp"
#include <cmath> // for fabs()
#include <utility> // for std::swap in orderVerticesCCWise()


namespace hms
{


/**
 * @brief It is often necessary to go from one end of an edge to the other.
 * This function makes that easier.
 * @param edgeVertices The table holding the vertex indices for each edge
 * @param edge Index of the edge in question
 * @param currentVertex Index of the vertex on one end
 * @return Index of the vertex on the other end of an edge
 */
Index otherVertex(
	const Array2Xi& edgeVertices,
	Index edge,
	Index currentVertex
){
	return aux::otherEntry(edgeVertices.col(edge), currentVertex);
}

namespace aux {
	Index otherEntry(
		const Ref<const Array2i>& col,
		Index entry
	){
		assert( hasEntry(col, entry) );

		return col[0] == entry ?
			col[1] :
			col[0];
	}
}


/**
 * @brief compares the vertex indices of two edges
 * @param mesh
 * @param edgeA
 * @param edgeB
 * @return the index of the shared (=connecting) vertex of two edges or else the 
 * error code signalCode::not_found
 */
Index sharedVertex( const Array2Xi& edgeVertices, Index edgeA, Index edgeB ){
	for ( auto a : edgeVertices.col(edgeA) )
		if ( hms::hasEntry(edgeVertices.col(edgeB), a) )
			return a;
	
	return signalCode::not_found;
}

Index sharedVertex( const Mesh& mesh, Index edgeA, Index edgeB ){
	return sharedVertex( mesh.edges().vertices(), edgeA, edgeB );
}


/**
 * @overload
 */
Index sharedVertex( const Edge& edgeA, const Edge& edgeB ){
	return sharedVertex( edgeA.mesh(), edgeA.index(), edgeB.index() );
}


/**
 * @brief
 * @param mesh The parent mesh
 * @return a std::vector with the same number of items as the mesh's
 * vertices. At each position, the edges connected to the corresponding
 * vertex are stored, as a std::vector<Index> (IdxVector).
 */
std::vector<IdxVector> edgesAtVertices(
	Index nVertices,
	const Array2Xi& edgeVertices
){
	std::vector<IdxVector> connectedEdges( nVertices );
	
	/* loop over edges, take edge index and store it for both vertices */
	for ( Index i=0; i < edgeVertices.cols(); ++i )
	{
		connectedEdges[edgeVertices(0,i)].push_back(i);
		connectedEdges[edgeVertices(1,i)].push_back(i);
	}
	
	return connectedEdges;
}

std::vector<IdxVector> edgesAtVertices( const Mesh& mesh )
{
	return edgesAtVertices(
		mesh.nVertices(),
		mesh.edges().vertices()
	);
}


/**
 * @brief Finds closed loops formed by an edge with edges of a higher index.
 * The higher-index-requirement allows for uniqueness of the found loops.
 * @param edgeVertices the vertex indices of all edges
 * @param edge the edge for which loops are to be found
 * @param edgesAtVertex output of edgesAtVertices(mesh)
 * @return a std::vector with up to two entries, which represent the loops
 * found. These loops are intended to be used as cells. In 2D space, an edge
 * can only be part of up to two cells, that is why the return object has up
 * to two entries. These are of type std::vector<Index> and hold the indices 
 * of edges which together form a closed loop. Only triangles and quadrilaterals
 * are considered, therefore these std::vector<Index> have length 3 or 4.
 */
std::vector<IdxVector> findLoops(
	const Array2Xi& edgeVertices,
	Index edge,
	const std::vector<IdxVector>& edgesAtVertex
){
	std::vector<IdxVector> loops;
	auto vertAtStart { edgeVertices(0,edge) };
	auto vertAtEnd   { edgeVertices(1,edge) };
	bool edgeAtStartUsed {false};
	
	for ( auto edgeAtStart : edgesAtVertex[vertAtStart] )
	{
	/* because the current edge is also among edgesAtVertex[vertAtStart],
	 * we need to guard against the current edge being found as the perfect
	 * loop. Additionally, edges for which cells were already found, i.e.
	 * all edges with a lower index than the current edge, must also be
	 * excluded. Otherwise, the same cell will be found 3-4 times.
	 */
		edgeAtStartUsed = false;
		if ( edgeAtStart <= edge ) continue;
		for ( auto edgeAtEnd : edgesAtVertex[vertAtEnd]   )
		{
			if ( edgeAtEnd <= edge ) continue;
			/* Stage 1: Find triangles.
			 * The edges connected to either end of the current edge need to 
			 * share a vertex to form a loop. */
			
			/* We start by getting the vertices at the other end of the 
			 * connected edges */
			auto a { otherVertex( edgeVertices, edgeAtStart, vertAtStart ) };
			auto b { otherVertex( edgeVertices, edgeAtEnd  , vertAtEnd   ) };
			
			if ( a == b ) // if these are the same, that's a triangle
			{
				loops.push_back( {edge, edgeAtEnd, edgeAtStart} );
				/* there can at maximum be two cells per edge */
				if (loops.size() == 2)
					return loops;
				/* guard against a triangle being found as a quadrilateral.
				 * breaking out of this loop goes to the next edgeAtStart */
				break;
			}
			/* Stage 2: Find quadrilaterals
			 * Two edges connected to the current edge must have another 
			 * connecting edge. */
			
			/* Now we retrieve all edges which are connected to A and B
			 * vertices and look for matches */
			for ( auto edgeA : edgesAtVertex[a] )
			{
				/* the edge closing the loop must also be unused thus far */
				if ( edgeA <= edge ) continue;
				for ( auto edgeB : edgesAtVertex[b] )
				{
					if ( edgeB <= edge ) continue;
					if ( edgeA == edgeB ) // quadrilateral found
					{
						loops.push_back(
							{ edge, edgeAtEnd, edgeA, edgeAtStart } );
						if (loops.size() == 2)
							return loops;
						/* there cannot be another loop formed with the same
						 * edge from the outermost loop, so we break all the 
						 * way through to the next one */
						edgeAtStartUsed = true;
						break;
					}
				}
				if (edgeAtStartUsed)
					break;
			}
			if (edgeAtStartUsed)
				break;
		}
	}
	
	return loops;
}

/**
 * @overload
 */
std::vector<IdxVector> findLoops(
	const Mesh& mesh,
	Index edge,
	const std::vector<IdxVector>& edgesAtVertex
){
	return findLoops( mesh.edges().vertices(), edge, edgesAtVertex );
}

/**
 * @overload
 */
std::vector<IdxVector>  findLoops(
	const Mesh& mesh,
	const Edge& edge,
	const std::vector<IdxVector>& edgesAtVertex
){
	return findLoops( mesh, edge.index(), edgesAtVertex );
}



/**
 * @brief Given that in a loop of edges, neighbouring edges are connected at a
 * shared vertex, the total vertex indices of the loop contain duplicates.
 * @param edgeVertices vertex indices of all edges
 * @param edgeLoop std::vector of the edge indices
 * @return a std::vector of the vertex indices forming the edge loop corners,
 * without duplicates, unordered.
 */
ArrayXi edgeLoopToVertices(
	const Array2Xi& edgeVertices,
	const Ref<const ArrayXi>& edgeLoop
){
	/* this slightly more involved heuristic approach is >20% faster than 
	 * collecting all vertices, sorting and removing duplicates, using 
	 * std::sort and std::unique */
	
	/* for a closed loop, there is the same number of edges and vertices */
	ArrayXi verts ( edgeLoop.size() );
	
	/* we know that thus far, both of the first edge's vertices are unique */
	verts.head(2) = edgeVertices.col( edgeLoop[0] );

	/* for triangles, only one more vertex is needed, which is shared by
	 * the other two vertices */
	if (edgeLoop.size() == 3)
	{
		verts[2] = sharedVertex( edgeVertices, edgeLoop[1], edgeLoop[2] );
		if (verts[2] == signalCode::not_found)
			throw hms::MeshError(
				"Loop of edges " + hms::toString( edgeLoop ) +
				" not closed!"
			);
	}
	else if (edgeLoop.size() == 4)
	{
		/* two vertices are known, the other ones are found by looking at
		 * the remaining edges and checking whether their vertices are
		 * included already */
		Index i {2};
		for ( char j {1}; j < edgeLoop.size(); ++j )
		{
			for ( Index vertex : edgeVertices.col( edgeLoop[j] ) )
			{
				if ( !(verts.head(i) == vertex).any() ){
					verts(i++) = vertex;
				}
				/* as soon as all vertices are found, return */
				if ( i == edgeLoop.size() ) return verts;
			}
		}
	}
	else // this should not happen
		throw hms::MeshError(
			"Loop of edges " + hms::toString( edgeLoop ) +
			" has unexpected size!"
		);
	
	return verts; /* to satisfy compiler */
}

/**
 * @overload
 */
ArrayXi edgeLoopToVertices(
	const Mesh& mesh,
	const Ref<const ArrayXi>& edgeLoop
){
	return edgeLoopToVertices( mesh.edges().vertices(), edgeLoop );
}


/**
 * @overload
 */
ArrayXi edgeLoopToVertices(const Mesh& mesh, const IdxVector& edgeLoop)
{
	return edgeLoopToVertices(
		mesh,
		Eigen::Map<const ArrayXi> ( edgeLoop.data(), edgeLoop.size() )
	);
}

ArrayXi edgeLoopToVertices( const Mesh& mesh, Index cell ){
	assert( !mesh.cell(cell).isGhost() );
	return edgeLoopToVertices(
		mesh,
		mesh.cell(cell).edges().head( mesh.cell(cell).nEdges() )
	);
}



/**
 * @brief takes a line AB and a point P, returns on which side of AB P lies.
 * All parameters are 2D cartesian vectors.
 * @param a beginning of line AB
 * @param b end of line AB
 * @param p point
 * @return an enum class object which represents the side of the line AB on
 * which the point p was found. Can be "left", "right" or "on". "on" is 
 * also returned, if the point is approximately on the line, within the set
 * relative epsilon value.
 */
Side whichSide(
	const Ref<const Vector2s>& a,
	const Ref<const Vector2s>& b,
	const Ref<const Vector2s>& p
){
	/* determinant method:
	 *         | 1  x_A  y_A |
	 * det M = | 1  x_B  y_B |
	 *         | 1  x_P  y_P |
	 * If det M < 0, P lies to the right, and if det M > 0 then P lies to
	 * the left of the line AB. If det M == 0, P lies on the line.
	 */
	MatrixNNs<3,3> m;
	m << 1, a.transpose(),
		 1, b.transpose(),
		 1, p.transpose();
	
	auto det { m.determinant() };
	if ( hms::approxEqual( det, 0 ) )
		return Side::on;
	return (det < 0 ? Side::right : Side::left );
}

/**
 * @overload
 */
Side whichSide( const Edge& e, const Ref<const Vector2s>& p ){
	return whichSide(
		e.mesh().vertex(e.vertices()[0]),
		e.mesh().vertex(e.vertices()[1]),
		p );
}

/**
 * @overload
 */
Side whichSide( 
	const Mesh& mesh,
	Index i,
	const Ref<const Vector2s>& p
){
	return whichSide(
		mesh.vertex(mesh.edges().vertices(i)[0]), // begin of edge
		mesh.vertex(mesh.edges().vertices(i)[1]), // end of edge
		p );
}



/**
 * @brief finds the index of the bottom left vertex among the passed-in vector. 
 * Lexicographic sorting, with y value first, then x value.
 * @return the column number of the bottom left vertex index within the
 * passed-in matrix
 */
Index findBottomLeft( const Ref<const Matrix2Xs>& vertices ){
	/* find lowest y value and its index */
	Index minInd;
	scalar minY { vertices.row(1).minCoeff( &minInd ) };
	scalar minX { vertices.col(minInd).x() };

	/* now we check whether there's a vertex with an approximately equal 
	 * y value but lower x value */
	for ( Index i{0}; i<vertices.cols(); ++i ){
		if (
			approxEqual( minY, vertices.col(i).y() ) &&
			vertices.col(i).x() < minX
		)
			minInd = i;
	}
	return minInd;
}

Index findBottomLeft( const Mesh& mesh, const ArrayXi& vertices ){
	return findBottomLeft( mesh.vertices()(hms::all, vertices) );
}


/**
 * @brief Orders vertex indices according to their cartesian coordinates as
 * either clockwise or counterclockwise. First index will always be that of the
 * bottom left vertex.
 * @param mesh parent mesh that holds the vertices
 * @param vertices vertex indices to be ordered
 * @param order can be either VertexOrdering::clockwise or 
 * VertexOrdering::counterclockwise. Default is counterclockwise.
 */
void orderVertices(
	const Matrix2Xs& vertexCoords,
	Ref<ArrayXi> vertices,
	VertexOrdering order
){
	/* First, find bottom left vertex (BL), then sort by (pseudo-)angle formed
	 * by line between BL and current vertex and x-axis
	 */
	std::swap(
		vertices[ 0 ],
		vertices[ findBottomLeft( vertexCoords(hms::all, vertices) ) ]
	);
	
	/* bubble sort vertices (std::sort with std::bind was veery slow?!) */
	bool sorted {false};
	while (!sorted)
	{
		sorted = true;
		for ( Index i {1} ; i < vertices.size()-1 ; ++i )
		{
			/* can't wait for C++20 to introduce "using enum"! */
			if ( pseudoangle(
				vertexCoords.col(vertices[i + (order==VertexOrdering::clockwise? 0 : 1 )])
				 - vertexCoords.col(vertices[0])
				 ) < pseudoangle(
				vertexCoords.col(vertices[i + (order==VertexOrdering::clockwise? 1 : 0 )])
				 - vertexCoords.col(vertices[0])
				 )
			)
			{
				sorted = false;
				std::swap( vertices[i], vertices[i+1] );
			}
		}
	}
}

void orderVertices(
	const Mesh& mesh,
	Ref<ArrayXi> vertices,
	VertexOrdering order
){
	orderVertices( mesh.vertices(), vertices, order );
}



/**
 * @brief Takes a closed loop of edges by their indices and an ordered list of
 * their vertices. Orders them according to the vertex list, i.e. if the 
 * vertices are ordered clockwise, then the returned list of edges is also
 * ordered clockwise.
 * @param edgeVertices
 * @param orderedVertices
 * @param edgeLoop
 */
void orderEdgeLoop(
	const Array2Xi& edgeVertices,
	const Ref<const ArrayXi>& orderedVertices,
	Ref<ArrayXi> edgeLoop
){
	assert( orderedVertices.size() == edgeLoop.size() );
	/* use ordered vertices in pairs, compare with given edge loop and sort
	 * (selection sort) */
	for ( Index i {0}; i < edgeLoop.size()-1; ++i )
	{
		for ( Index j {i+1}; j < edgeLoop.size(); ++j )
		{
			if ( hms::sameEntries(
				edgeVertices.col( edgeLoop[j] ),
				orderedVertices.segment(i,2)
			) )
			{
				std::swap( edgeLoop[i], edgeLoop[j] );
				break;
			}
		}
	}
}

/**
 * @overload
 */
void orderEdgeLoop(
	const Mesh& mesh,
	const Ref<const ArrayXi>& orderedVertices,
	Ref<ArrayXi> edgeLoop
){
	orderEdgeLoop( mesh.edges().vertices(), orderedVertices, edgeLoop );
}

/**
 * @overload
 */
void orderEdgeLoop(
	const Mesh& mesh,
	const Ref<const ArrayXi>& orderedVertices,
	IdxVector& edgeLoop
){
	orderEdgeLoop(
		mesh,
		orderedVertices,
		Eigen::Map<ArrayXi>( edgeLoop.data(), edgeLoop.size() )
	);
}
	


IdxVector connectingEdges(
	const Mesh& mesh,
	const IdxVector& vertices,
	const std::vector<IdxVector>& edgesAtVertex
){
	using hms::hasEntry;
	/* Go through the edges connected to each specified vertex and see whether
	 * the other vertex of that edge is also part of the specified vertex list.
	 * This should find each edge twice, so guard against duplicates */
	IdxVector edges;
	
	for ( auto vertex : vertices )
	{
		for ( auto edge : edgesAtVertex[vertex] )
			if (
				hasEntry(vertices, otherVertex(mesh.edge(edge), vertex) ) &&
				!hasEntry(edges, edge) ) // no duplicates
				edges.push_back(edge);
	}
	return edges;
}


AlignedBox boundingBox( const Ref<const Matrix2Xs>& polygon )
{
	/* bottom left and top right corner */
	return { polygon.rowwise().minCoeff(), polygon.rowwise().maxCoeff() };
}


AlignedBox boundingBox( const std::vector<Vector2s>& polygon ){
	return boundingBox(
		Eigen::Map<const Matrix2Xs>( &(polygon[0][0]), 2, polygon.size() )
	);
}





/**
 * @brief Tests whether a point lies within a polygon
 * 
 * @param polygon the corners of the polygon, in order (i.e. following along the
 * edges of the polygon)
 * @param point 
 * @return true if the point lies within the polygon
 * @return false if the point lies outside the polygon
 * 
 *  
 * follows efficient winding number algorithm,
 * source: http://geomalgorithms.com/a03-_inclusion.html
 * Therefore, returns true for self-overlapping regions, and false for cutouts.
 */
bool inPolygon(
	const Ref<const Matrix2Xs>& polygon,
	const Ref<const Vector2s>& point
){
	int wn {0}; // winding number counter
	
	for ( Index j {0}; j < polygon.cols() ; ++j )
	{
		const Ref<const Matrix2Xs>::ConstColXpr nextCorner {
			polygon.col( (j+1) % polygon.cols())
		};

		if ( polygon(ydim,j) <= point[ydim] ){
			if ( nextCorner[ydim] > point[ydim] ) // potential upward crossing
				if ( whichSide( polygon.col(j), nextCorner, point ) == Side::left )
					++wn;
		} else {
			if ( nextCorner[ydim] <= point[ydim] )// potential downward crossing
				if ( whichSide( polygon.col(j), nextCorner, point ) == Side::right)
					--wn;
		}
	}
	return wn != 0; // when the winding number is 0, the point is outside
}

/**
 * @overload
 */
bool inPolygon(
	const std::vector<Vector2s>& polygon,
	const Ref<const Vector2s>& point
){
	return inPolygon(
		Eigen::Map<const Matrix2Xs>(&(polygon[0][0]), 2, polygon.size() ), point
	);
}



IdxVector verticesInPolygon(
	const Mesh& mesh,
	const std::vector<Vector2s>& polygon
){
	IdxVector verticesInside;
	
	/* no vertex can be inside a line */
	if (polygon.size() < 3)
		return {};
	
	AlignedBox boundingBox { hms::boundingBox(polygon) };
	
	for ( Index i {0}; i < mesh.nVertices() ; ++i )
	{
		/* exclude those points which aren't inside the polygon's bounding box*/
		if ( !boundingBox.contains(mesh.vertex(i)) )
			continue;
		if ( inPolygon( polygon, mesh.vertex(i) ) )
			verticesInside.push_back(i);
	}
	return verticesInside;
}

IdxVector edgesInPolygon( 
	const Mesh& mesh,
	const std::vector<Vector2s>& polygon,
	const std::vector<IdxVector>& edgesAtVertex
){
	return connectingEdges(
		mesh, verticesInPolygon(mesh, polygon),
		edgesAtVertex.empty() ?
		edgesAtVertices(mesh) :
		edgesAtVertex
	);
}

IdxVector cellsInPolygon(
	const Mesh& mesh,
	const std::vector<Vector2s>& polygon
){
	IdxVector cellsInside;
	
	/* no vertex can be inside a line */
	if (polygon.size() < 3)
		return {};
	
	AlignedBox boundingBox { hms::boundingBox(polygon) };
	
	/* this would NOT work for structured meshes with simplified indexing! */
	for ( Index i = 0; i<mesh.nCells(); ++i )
	{
		if ( !boundingBox.contains( mesh.cell(i).centroid() ) )
			continue;
		if ( inPolygon( polygon, mesh.cell(i).centroid() ) )
			cellsInside.push_back(i);
	}
	return cellsInside;
}



IdxVector verticesInBox(
	const Mesh& mesh,
	const Vector2s& a,
	const Vector2s& b
){
	IdxVector verticesInside;
	
	AlignedBox boundingBox { hms::boundingBox( {a, b} ) };
	
	for ( Index i {0} ; i < mesh.nVertices() ; ++i )
	{
		if ( boundingBox.contains( mesh.vertex(i) ) )
			verticesInside.push_back(i);
	}
	return verticesInside;
}

IdxVector edgesInBox(
	const Mesh& mesh,
	const Vector2s& a,
	const Vector2s& b,
	const std::vector<IdxVector>& edgesAtVertex
){
	return connectingEdges(
		mesh,
		verticesInBox(mesh, a, b),
		edgesAtVertex.empty() ?
		edgesAtVertices(mesh) :
		edgesAtVertex
	);
}


IdxVector cellsInBox(
	const Mesh& mesh,
	const Vector2s& a,
	const Vector2s& b
){
	IdxVector cellsInside;
	
	auto boundingBox { hms::boundingBox( std::vector{a, b} ) };
	
	/* this would NOT work for structured meshes with simplified indexing! */
	for ( Index i = 0; i<mesh.nCells(); ++i )
		if ( boundingBox.contains( mesh.cell(i).centroid() ) )
			cellsInside.push_back(i);
	
	return cellsInside;
}



IdxVector verticesToEdgeLine(
	const Mesh& mesh,
	const IdxVector& vertices, 
	const std::vector<IdxVector>& edgesAtVertex
){
	IdxVector edgeLine;
	
	/* go through specified vertices and check out each edge connected to
	 * each vertex. If both of the edge's vertices are in the specified list,
	 * then it's part of the line. */
	for ( auto vertex : vertices )
	{
		for ( auto edge : edgesAtVertex[vertex] )
		{
			if (
				hms::hasEntry(
					vertices,
					otherVertex(mesh.edge(edge), vertex) ))
			{
			/* To prevent duplicates, check last and second last entry */
				if (edgeLine.size()>=1 && edgeLine.back() == edge )
					continue;
				if (edgeLine.size()>=2 && edgeLine[edgeLine.size()-2] == edge )
					continue;
				
				edgeLine.push_back(edge);
			}
		}
	}
	
	return edgeLine;
}


template<typename MeshType>
std::optional<Index> containingCell(
	const MeshBase<MeshType>& meshBase,
	Array2s coords
){
	const MeshType& mesh { meshBase.derived() };
	[[maybe_unused]] static constexpr Index ng {1};

	if constexpr (std::is_same_v<MeshType, UniMesh>){
		if (
			(coords >= mesh.origin().array() - (ng * mesh.diff()) ).all() &&
			(coords <= mesh.origin().array() +
				(mesh.size() + ng).template cast<scalar>() * mesh.diff()
			).all()
		){
			Array2i ij { ( (coords - mesh.origin().array() ) / mesh.diff() )
				.template cast<Index>()
			};
			return toSingleIndex(mesh, ij[0], ij[1]);
		}
	} else if constexpr (std::is_same_v<MeshType, RectMesh>){
		auto
			beg_x { mesh.xCoords().cbegin() },
			beg_y { mesh.yCoords().cbegin() },
			end_x { mesh.xCoords().cend() },
			end_y { mesh.yCoords().cend() },
			it_x  { std::upper_bound( beg_x, end_x, coords.x() ) },
			it_y  { std::upper_bound( beg_y, end_y, coords.y() ) };
		if (
				it_x > beg_x && it_x < end_x
			&& it_y > beg_y && it_y < end_y
		)
			return toSingleIndex( mesh,
				it_x - (beg_x + ng + 1),
				it_y - (beg_y + ng + 1)
			);
	} else {
		throw std::runtime_error("containingCell not implemented for non-rectangular meshes!");
	}
	return std::nullopt;
}


template<typename MeshType>
IdxVector cellsAlongLine( const MeshBase<MeshType>& mesh,
	Array2s lineBeg,
	Array2s lineEnd
){
	mesh.logger().trace(
		"Tracing cells between coords {} and {}...",
		lineBeg.transpose(), lineEnd.transpose()
	);
	Matrix2s line, verts;
	line << lineBeg.matrix(), lineEnd.matrix();

	std::optional<Index>
		cellBeg { containingCell(mesh, std::move(lineBeg) ) },
		cellEnd { containingCell(mesh, std::move(lineEnd) ) };

	if ( !cellBeg.has_value() || !cellEnd.has_value() ){
		throw std::runtime_error("Line start or end outside of mesh!");
	} else {
		mesh.logger().trace("Start and end cell indices are {} and {}.",
			cellBeg.value(), cellEnd.value()
		);
	}

	using CellType = typename MeshType::CellType;
	using EdgeType = typename MeshType::EdgeType;

	Index currentIndex { cellBeg.value() }, prevIndex, nextIndex {currentIndex};

	IdxVector cells;
	bool deadEnd {false};

	/* go along the line, check which edges it intersects 
	 * and use the cell on the other side of the edge as the new current cell,
	 * until the last cell is reached */
	while( !deadEnd ){
		/* prevent infinity loop by going to the next cell, then back and forth */
		prevIndex = currentIndex;
		currentIndex = nextIndex;
		cells.push_back( currentIndex );
		mesh.logger().trace( "Adding cell #{} to line.", currentIndex );

		if ( currentIndex == cellEnd.value() )
			break;

		CellType cell { mesh.cell(currentIndex) };

		if ( cell.isGhost() )
			break;

		/* each edge of the cell is checked for an intersection,
		 * and if one is found, go to the cell on the other side */
		deadEnd = true;
		for ( Index edgeIndex : cell.edges() ){
			EdgeType edge { mesh.edge(edgeIndex) };
			verts <<
				mesh.vertex( edge.vertices()[0] ),
				mesh.vertex( edge.vertices()[1] );

			if ( intersects(line, verts) ){
				const Ref<const Array2i>& adjCells { edge.cells() };
				nextIndex = adjCells[0] == currentIndex ? adjCells[1] : adjCells[0];
				if ( nextIndex != prevIndex ){
					deadEnd = false;
					break;
				}
			}
		}
	}
	return cells;
}


#define INSTANTIATE_FOR_ALL_MESHTYPES(MeshType) \
template std::optional<Index> containingCell(const MeshBase<MeshType>&, Array2s); \
template IdxVector cellsAlongLine( const MeshBase<MeshType>&, Array2s, Array2s );

#include "meshInstantiator.hpp"

} // end namespace hms